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ABSTRACT 

Several observations of transition discs show lopsided dust-distributions. A potential 
explanation is the formation of a large-scale vortex acting as a dust-trap at the edge of 
a gap opened by a giant planet. Numerical models of gap-edge vortices have so far em¬ 
ployed locally isothermal discs in which the temperature profile is held fixed, but the 
theory of this vortex-forming or ‘Rossby wave’ instability was originally developed for 
adiabatic discs. We generalize the study of planetary gap stability to non-isothermal 
discs using customized numerical simulations of disc-planet systems where the planet 
opens an unstable gap. We include in the energy equation a simple cooling function 
with cooling timescale tc = where flk is the Keplerian frequency, and examine 

the effect of /3 on the stability of gap edges and vortex lifetimes. We find increas¬ 
ing P lowers the growth rate of non-axisymmetric perturbations, and the dominant 
azimuthal wavenumber m decreases. We find a quasi-steady state consisting of one 
large-scale, over-dense vortex circulating the outer gap edge, typically lasting O(IO^) 
orbits. We find vortex lifetimes generally increase with the cooling timescale tc up 
to an optimal value of tc 10 orbits, beyond which vortex lifetimes decrease. This 
non-monotonic dependence is qualitatively consistent with recent studies using strictly 
isothermal discs that vary the disc aspect ratio. The lifetime and observability of gap- 
edge vortices in protoplanetary discs is therefore dependent on disc thermodynamics. 

Key words: accretion, accretion discs, protoplanetary discs, hydrodynamics, insta¬ 
bilities, planet-disc interactions, methods: numerical 


1 INTRODUCTION 


The interaction between planets and protoplanetary discs 
plays an important role in the theory of planet formation 
and disc evolution. Disc-planet interaction may lead to the 
orbital migration of protopl anets and modify the stru cture 
of protoplanetary discs (see iBaruteau &: Massetll2013l . for a 
recent review). 

A sufficiently massive planet 
p rotoplanetary 


open 




_ gaseous protoplanetary disc lPapaIoizpu_J^Lin| 

Bg;den et ajJ[]j99|: ICrida. Morbidelli fc Massetll2006l: 


Fung. Shi &i Chiansll2014 ). while low mass pl anets may also 
open gaps if the disc viscosity i s small enough (IlT et al. 200d; 
iDong. Rafikov Stond l201ll : iDuffell fc MacFadverJ 20131) . 
Support for such disc-planet interaction have begun to 
emerge in observa tions of circumstella r discs that reveal 
annular gaps (e.g. lOuanz et al. I l2013bl: iDebes et al.l l2013l: 
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lOsorio et al. I B). with possible evidence of companion s 
within them (e.g. Touanz et al.ll201^ : iReggiani et al.l[^014l ) . 

A recent theoretical development in the study of plan¬ 
etary gaps is their stability. When the disc viscosity is low 
and/or the planet mass is large, the presence of potential 
vorticity (PV, the ratio of vorticity to surface density) ex¬ 
trema can render planetary gaps dynamically unstable due 
to wh at is now referred to as the ‘Rossby wave instability’ 


(RWI. 1 Lovelace et al.l If 9991: 

Li et al.l 

200c 

). This eventually 

leads to vortex formation ( 

Li et al.l 

2001 

: iKoller. Li & LirJ 

I2OO3I: iLi et al.l l2005l: Ide Val-Borro et al.l 

2007). which can 


significantly affect orbital migration of the plane t JOu e t al.l 

I2OO7I : iLi et al.ll2009l : IYu et ah I2OI0I : iLin fc PaDaloizoull201Cll ). 


Vortex formation at gap edges may also have observ¬ 
able consequences. Because disc vortices represent pres- 


sure 

jBari 

maxima, the^ 

ge & Sommerial 

r are able to collect dust particles 
Il995l:llnaba & Bargell2006l:lLvra & LirJ 

|2013t 

). Dust-trapping at gap-edge vortices have thus been 


suggested to explain asymmetric d ust distributions ob- 
served in several transition discs (e.g. ICasassus et al. I l2013l : 
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vandei_^axd_et_^ ]|2013l: IseUa^^t_al I gOjj : iFukagawa et all 


2 OI 3 I : Perez et al. 2014 : Pinilla et al. l2015h . 


However, studies of Rossby vortices at planetary gap- 
edges have adopted locally isothermal discs, where the 
disc temperature is a fixed function of position only (e.g . 


Lvra^t^l. 20091 : iLin fc Papaloizoul I 2 OIII : IZhu et 


p osition only (e.g . 
IZhu et ^I20l4 


Fu et al. 2 OI 4 I I. On the other hand, the the ory of th e RWI 


was in fact developed for adiabatic discs dLi et al.l [200011 ■ 
which permits heating. In adiabatic discs, the relevant quan¬ 
tity for stability becomes a ge neralization of the P V that 
accounts for entropy variations dLovelace et al.iri999h . 

Gap-opening is associated with planet-induced spiral 
shocks. In an isothermal disc, PV- generation across these 
isothermal shocks leads to the RW I (|K(hler^i_J^Jhn||200^; 
Li et al.l I 2 OO 5 I : Ide Val-Borro et akl I 2 OO 7 I : Lin fc Papaloizoul 


2 OIOII . However, if cooling is inefficient and the shock is non- 


isothermal, then shock-heating may affect gap stability, since 
the relevant quantity is an entropy-modified PV (described 
below), and there is entropy-generation across the shocks. 

For example, previous linear simulations of the RWI 
found that increasing t he sound-speed favours instability 
dLi et al.ll200Cll : lLirJl2013l L In the context of planetary gaps, 
however, the increased temperature may also act to stabilize 
the disc by making gap-opening more difficult. It is there¬ 
fore of theoretical interest to clarify the effect of heating and 
cooling on the stability of planetary gaps. 

In this work, we extend the study of planetary gap sta¬ 
bility against vortex formation to non-isothermal discs. We 
include in the fluid energy equation an one-parameter cool¬ 
ing prescription that allows us to probe disc thermodynam¬ 
ics ranging from nearly isothermal to nearly adiabatic. 

This paper is organized as follows. In ^Uwe describe the 
equations governing the disc-planet system and initial con¬ 
ditions. Our numerical approach, including diagnostic mea¬ 
sures, are given in ^ We present results from two sets of 
numerical experiments. In )2]we use disc-planet interaction 
to set up discs with gaps, but study their stability without 
further influence from the planet. We then perform long¬ 
term disc-planet simulations to examine the lifetime of gap- 
edge vortices in an as a function of the imposed cooling 
rate. We conclude and summarize in aSJwith a discussion of 
important caveats. 


2 DISC-PLANET MODELS 

The system is a two-dimensional (2D) non-self-gravitating 
gas disc orbiting a central star of mass M«. We adopt 
cylindrical co-ordinates (r, (j), z) centred on the star. The 
frame is non-rotating. Computational units are such that 
G = M, — TZ = — 1 where G is the gravitational con¬ 

stant, TZ is the gas constant and p is the mean molecular 
weight. 

The disc evolution is governed by the standard fluid 
equations 

BT 

^ + V • (Er;) = 0, (1) 

^+„.Vr; = -ivp-V-I> + /., (2) 

Oc 

— +\/-{ev) = -pV-v + H-C, (3) 


where E is the surface density, v = {vr, v^) the fluid velocity, 
p is the vertically-integrated pressure, e = p /(7 — 1 ) is the 
energy per unit area and the adiabatic index 7 = 1.4 is 
assumed constant. 

The total potential <I> includes the stellar potential, 
planet potential (described below) and indirect potentials to 
account for the non-inertial reference frame. In the momen¬ 
tum equations, represent viscous forces, which includes 
artificial bulk viscosity to handle shocks, and a Navier- 
Stokes viscosity whose magnitude is characterized by a con¬ 
stant kinematic viscosity parameter v. However, we will be 
considering effectively inviscid discs by adopting small val¬ 
ues of V. 


2.1 Heating and cooling 

In the energy equation, the heating term H is defined as 

n = Q+-Qt^, ( 4 ) 


where Q'^ represents viscous heating (from both physical 
and artificial viscosity) and subscript i denotes evaluation 
at t = 0. The cooling term C is defined as 



where tc = is the cooling time, GM/r^ is 

the Keplerian frequency and /3 is an input parameter. This 
cooling prescription allows one to explore the full range of 
thermodynamic response of the disc in a systematic way: 
/3 1 is a locally isothermal disc while /I 3 > 1 is an adiabatic 

disc. 

Note that the energy source terms have been chosen 
to be absent at t = 0 , allowing the disc to be initialized 
close to steady state. The C function attempts to restore the 
initial energy density (and therefore temperature) profile. In 
practice, this is a cooling term at the gap edge because disc- 
planet interaction leads to heating. 


2.2 Disc model and initial condition 

The disc occupies r £ [nn,rout] and cf) £ [0, 27r]. The initial 
disc is axisymmetric with surface density profile 


EW = E,.,(^) 



rin 

r -I- Hi{rin) 


( 6 ) 


where the power-law index s = 2, H{r) = CisoH/; defines 
the disc scale-height where Ciso = \fpfL is the isothermal 
sound-speed. The disc aspect ratio is defined aa h = H/r 
and initially h = 0.05. For a non-self-gravitating disc, the 
surface density scale Eref is arbitrary. 

The initial azimuthal velocity v^i is set by centrifugal 
balance with pressure forces and stellar gravity. For a thin 
disc, ~ rQk- The initial radial velocity is Vr = Sv/r, 
where u = i>rf„Hj;(rin), and we adopt 0 = 10 “®, so that 
\vr/vtt,\ 1 and the initial flow is effectively only in the 
azimuthal direction. With this value of physical viscosity, 
the only source of heating is through compression, shock¬ 
heating (via artificial viscosity) and the C function when 
e/E < ei/Ei. 


© 0000 RAS, MNRAS 000, 000-000 




















































Gaps in non-isothermal discs 3 


2.3 Planet potential 


The planet potential is given by 


= - 


GM„ 


7|T--rp|2 + e2 


(7) 


where Mp is the planet mass and we fix g = Mp/Mt = 
10“® throughout this work. This corresponds to a Jupiter- 
mass planet if M* = Mq. The planet’s position in the disc 
Vp = {rp,(f>p) and tp = 0.5rh is a softening length with 
rh ~ being the Hill radius. The planet is held on 

a fixed circular orbit with Tp = lOrin and 4>p = Q.k{rp)t. This 
also defines the time unit Po = 2n/Qk{rp) used to describe 
results. 


3 NUMERICAL EXPERIMENTS 


The disc-planet system is ev olved using the FARGO-ADSG code 
dSaruteau fc Massefjl200 8^ This is a modified version of 
the original FARGO code ( Masset|[2000li to include the energy 
equation. The code e mploys a finite-difference scheme sim¬ 
ilar to the ZEUS code Istone &: Normanlll992li . but with a 
modified azimuthal transport algorithm to circumvent the 
time-step restriction set by the fast rotation speed at the 
inner disc boundary. The disc is divided into {Nr, Nrj,) zones 
in the radial and azimuthal directions, respectively. The grid 
spacing is logarithmic in radius and uniform in azimuth. 


3.1 Cooling prescription 

In this work we only vary one control parameter: the cooling 
time. The cooling parameter /3 is chosen indirectly through 
the parameter P such that 


tr{rp + Xs) = ^(rp -I- Is) = /3tiib(rp + is), (8) 


where is is the distance from the planet to its gap edge, and 
tiib is the time interval between successive encounters of a 
fluid element at the gap edge and the planet’s azimuth. That 
is, we measure the cooling time in units of the time interval 
between encounters of a fluid element at the gap edge and 
the planet-induced shock. 

Assuming Keplerian orbital frequencies and is ^ ip 
gives tiib ~ 47rrp/(3flfcpis), where Q,kp = flk{rp). Therefore 


P 




(9) 


where is -C ip was used again. We use is = 2rh in Eq.[9l For 
a planet mass with q = 10~®, Eq. [9]then gives /3 ~ 23.9/3. 
In terms of planetary orbital periods, this is 


tcir) 


£ 

2-k 


/ \ 3/2 

(u) 



( 10 ) 


3.2 Diagnostic measures 

3.2.1 Generalised potential vorticity 

The generalised potential vorticity is defined as 

— J1 _ X 
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( 11 ) 


where = r- '^dr{r^Q.^) is the square of the epicyclic fre¬ 
quency, Q. = v^jr is the angular speed, and S = p/TP is the 
entropy. The first factor is the usual potential vorticity (PV, 
or vortensity). 

The generalised PV appears in the description of 
the linear stab i lity of radially-struc tured adiabatic discs 
llLovelace et al.l Il999l : iLi et al.l I2000I1 . where the authors 
show an extremum in fj may lead to a dynamical instability, 
the RWI. In a barotropic disc where p = p{T), the entropy 
factor is absent and the important quantity is the PV. 


3.2.2 Fourier modes 

The RWI is characterized by exponentially growing pertur¬ 
bations. Though in this paper we do not consider a formal 
linear instability calculation, modal analysis will be useful 
to analyse the growth of perturbations with different az¬ 
imuthal wavenumbers, which is associated with the number 
of vortices initially formed by the RWI. 

The Fourier transform of the time-dependent surface 
density is 

d(l) (12) 

where m is the azimuthal wave number. We define the 
growth rate a of the component of the surface density 
through 

(13) 

where {| • |)r denotes the average of the absolute value over 
a radial region of interest. By using Eq.[T3]the growth rates 
of the unstable modes can be found from successive spatial 
Fourier transforms over an appropriate period of time. 


3.2.3 Rossby number 


The Rossby number 


Ro = 


z-Vxu — (z-Vx v)^ 


(14) 


is a dimensionless measure of relative vorticity. Here {■)^ de¬ 
notes an azimuthal average. Values of Ro < 0 correspond to 
anti-cyclonic rotation with respect to the background shear 
and thus can be used to identify vortices and quantify its 
intensity. 


4 GROWTH OF NON-AXISYMMETRIC 
MODES WITHOUT THE INFLUENCE OF 
THE PLANET 

In this section, the planet is introduced a.t t — 20Po and its 
potential is switched on over IOPq. At t = 30Po we switch 
off the planet potential and azimuthally average the surface 
density, energy and velocity fields. (At this point the planet 
has carved a partial gap and the RWI has not yet occurred.) 
Effectively, we initialise the disc with a gap profile. We then 
perturb the surface density in the outer disc (r > ip) and 
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continue to evolve the disc. We impose sinusoidal pertur¬ 
bations with azimuthal wavenumbers m G [1,10] and ran¬ 
dom amplitudes within ±0.01 times the local surface den¬ 
sity. This procedure allows us to analyse the growth of non- 
axisymmetric modes associated with the gap, but without 
complications from non-axisymmetry arising directly from 
disc-planet interaction (i.e. planet-induced wakes). 

Note that these ‘planet-off’ simulations are not linear 
stability calculations because the cooling term in our energy 
equation restores the initial temperature profile correspond¬ 
ing to constant H/r — 0.05, rather than the heated gap edge. 
However, we will examine a nearly adiabatic simulation in 
114.3.11 which is closer to a proper linear problem. 

Simulations here employ a resolution of {Nr,N^) = 
(1024, 2048) with open boundaries at r = nn and rout = 
25rin. We compare cases with (3 = 0.1,1,10 corresponding 
to fast, moderately, and slowly cooled discs. 


4.1 Gap structure 

We first compare the gap structures formed by planet- 
disc interaction as a function of the cooling time. The 
azimuthally-averaged gap profiles are shown in Fig. [T] for 
different values of (3. Gaps formed with lower j3 (faster cool¬ 
ing) are deeper with steeper gradients at the gap edges. 
Faster cooling rates also increase (decrease) the surface den¬ 
sity maxima (minima). However, a clean gap does not form 
in this short time period. 

Increasing {3 leads to higher disc aspect ratios h = H/r, 
i.e. higher temperatures. Heating mostly occurs at the gap 
edges due to planet-induced spiral shocks. Increasing the 
cooling timescale implies that this heat is retained in the 
disc. In the inviscid li mit the gap opening condition i s 
rj, ^ H or q > 3/i® dCrida. Morbidelli fc Massed 1200^ . 
which indicates that for hotter discs (higher h), it becomes 
more difhcult for a planet of fixed q to open a gap. This 
explains the shallower gaps in surface density when j3 is in¬ 
creased. 

The important consequence of a heated gap edge is that 
the generalised vortensity profiles, fj, becomes smoother with 
increasing cooling times, with the extrema becoming less 
pronounced. Previous locally isothermal disc-plane t simula¬ 
tions show the RWI associat ed with PV minima (iLi et al.l 
l2005l : lLin fc PaDaloizoull2010l l. We can therefore expect the 
RWI to be associated with minima in the generalised vorten¬ 
sity (corresponding to local surface density maxima) in the 
non-isothermal case. Because the extrema are less sharp, the 
RWI is expected to be weaker and the gap to be more stable 
with longer cooling times. 

However, we remark that the change in the gap struc¬ 
ture becomes less significant at long cooling times, as Fig.[T] 
shows that the prohles with /3 = 1 and /? = 10 are similar. 
This implies that the effect of cooling timescale on the RWI 
through the set up of the gap profile, becomes less important 
for large (J. 


4.2 Axisymmetric stability 

The initial planet-disc interaction form bumps and grooves 
in the gap profiles which can potentially be unstable due 



-10 -50 5 10 

(r-rj/h 


Figure 1. Azimuthally averaged gap profiles at t = SOPq for the 
initial partial gap opened before instability emerges for fast (solid, 
p = 0.1), moderate (dashed, /3 = 1), and slow cooling (dashed- 
dot, [3 = 10). The relative surface density perturbation (top), 
disc aspect ratio (middle) and generalised vortensity perturbation 
(bottom) are shown. 

to axisymmetric instabilities. The generalised local axisym¬ 
metric stability condition is the Solberg-Hoiland criterion, 

± > 0 (15) 

where 


1 dP I 

r 1OE 

1 

dP\ 

E dr ' 

VE dr 

7P 

dr ) 


is the square of the Brunt-Vaisala frequency. At the outer 
gap edge r = rp ± 2.5rh, where the RWI is excited (see 
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Gaps in non-isothermal discs 5 


Table 1. Dominant mode and growth rates for /3 = 0.1,1.0,10.0 
(fast, moderate, and slow cooling) values during ‘planet-off’ sim¬ 
ulations 



0 

II 


p = 1.0 


p = 10.0 

m 

10^a-/Q(rp) 

m 

b 

(N 

0 

m 

10^a-/Sl(rp) 

6 

7.3 

3 

2.0 

1 

1.1 

7 

7.8 

4 

2.2 

2 

1.6 

8 

7.9 

5 

2.3 

3 

1.7 

9 

7.9 

6 

1.6 

4 

1.2 

10 

6.8 

7 

1.1 

5 

0.1 


,S=10.0 



t/Po 


;8 = 0.1, t = 50Po ;8=1.0, t = 67Po ^=10.0, t = 75Po 



2.5 3.0 3.5 2.5 3.0 3.5 2.5 3.0 3.5 

(r-g/r„ (r-g/r„ (r-g/r„ 


Figure 3. Generalised vortensity perturbation (relative to t = 0) 
for cases of /3 = 0.1,1,10 (left,middle,right) during the growth of 
non-axisymmetric modes. The planet potential has been switched 
off. The number of vortices decrease as /3 increases. Note that 
snapshots are taken later for increasing /3 because it takes longer 
for the vortices to grow and become visible with increasing cooling 
time. 


Figure 2. Evolution of azimuthal Fourier modes of disc surface 
density, non-dimensioqnalized by the initial axisymmetric com¬ 
ponent Uo(t = 0 ), for the ‘planet-off’ simulation with /3 = 10 . 
Colours correspond to different m values. The m = 3 component 
is the fastest growing mode during linear growth with a growth 
rate of 7 = 0.017f2(rp). 

below), we find -P reaches local minimum with a value 
~ 0.45 n^(rp) for all /3. The Brunt-Vaisala frequency at the 
outer gap edge is V ~ 0.1 fi(rp), decreasing marginally with 
longer cooling rate. The Solberg-Hoiland criteria is similarly 
satisfied for the entire 2D disc throughout the simulations. 
Thus for all values of /3 the planet-induced gaps are stable 
to axisymmertic perturbations. 

4.3 Non-axisymmetric instability 

We now examine the evolution of the gap for t > 30Po, 
with the planet potential switched off, but with an added 
surface density perturbation. For all three cooling times 
P = 0.1, 1.0, 10, we observe exponential growth of non- 
axisymmetric structures. An example is shown in Fig. for 
P = 10. We characterize these modes with an azimuthal 
wavenumber m and growth rate a{m) as defined by Eq.[T2]— 
Mode amplitudes were averaged over r — rp £ [2,5]r;i. 
Table[T]lists the growth rates measured during linear growth 
for 5 values of m centred around that with maximum growth 
rate. 

Table [T] show that as P is increased from 0.1 —^ 10 the 
dominant azimuthal Fourier mode decreases from m = 9 —^ 
3 and the respective growth rate decreases from 7 /D(rp) = 


0.079 —^ 0.017. However, despite two orders of magnitude 
increase in the cooling time, the instability remains dynam¬ 
ical with characteristic growth time < lOPo- Snapshots of 
the instability in for the different P are shown in Fig[3l 

These ‘planet-off’ simulations show that gap edges be¬ 
come more stable with longer cooling times. This is expected 
because larger P results in hotter gap profiles at t = 30Po 
with less pronounced generalised vortensity minima. Stabi¬ 
lization with increased cooling time is therefore due to a 
smoother basic state for the instability, as it is more diffi¬ 
cult for the planet to open a gap if the disc is allowed to 
heat up. 

For completeness we also simulated a locally isothermal 
disc with h = 0.05 where the sound-speed is kept strictly 
equal to its initial value. This simulation yield a most un¬ 
stable growth rate a/^{rp) ~ 0.05 at m = 5, compared with 
a value of 0.085 at m = 6 for a corresponding simulation 
that includes the energy equation but with rapid cooling 
P = 0.01. However, the vortex evolution is similar. 


4-3.1 Nearly adiabatic discs 

The above ‘planet-off’ simulations are not formally linear 
stability calculations, because the cooling time is compara¬ 
ble or shorter than the instability growth time, tc < 7~^. 
Thus the disc cools back to its initial temperature corre¬ 
sponding to h = 0.05 before or during the instability growth, 
so we do not have a steady basic state to formulate a stan¬ 
dard linear stability problem. 

In order to capture the effect of a heated gap edge, we 
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ran a simulation with $ = 100, corresponding to an almost 
adiabatic disc. In this simulation the cooling rate is slow 
enough that the gap temperature profile (e.g. middle panel 
of Fig. [T]) changes only marginally over the instability growth 
timescale. 

We find very similar gap profiles and mode growth rates 
for P = 100 as with ,0 = 10. At f = 30Po, the disc only 
heats up to values h ~ 0.06 in the nearly adiabatic case. 
This is close to the original temperature of h — 0.05, so 
linear growth r ates are not expected to change significantly 
dLi et al.ll2Oo'0'l . 

According to lLi et ^ (l2000l l. increasing h increases lin¬ 
ear growth rates of the RWI because it is pressure-driven. 
However, in the case of disc-planet interaction, increasing h 
has a stabilizing effect through the setting up the gap profile 
because it results in smoother gap edges. The fact that we 
observe smaller growth rates as h is increased indicates that 
for planetary gaps, the importance of h on the linear RWI 
is through setting up the gap profile, i.e. basic state for the 
instability (as opposed to the linear response). 

4.4 Long term evolution 

We also extended these ‘planet-off’ simulations into the non¬ 
linear regime. After the linear growth phase of the vortices, 
vortex merging takes hold on timescales of up to 150Po, until 
there is one vortex left. We find the vortex merging time is 
dependent on the growth rates of the modes and saturation 
timescales, with the slowest growing modes in (3 = 10 taking 
the longest to merge. 

FigSlshows evolution of the m = 1 surface density com¬ 
ponent, which represents the amplitude of the post-merger 
single vortex. For completeness we also ran intermediate 
cases with P — 0.5 and 5.0. The amplitude of the initially 
formed vortex was found to decrease with increased cooling 
rate. The vortices simply decay on a timescale of 0(10®) 
orbits with faster decay for stronger vortices (which are ob¬ 
tained with faster cooling rates). 

This decay is probably due to numerical viscosity. Dur¬ 
ing the slow decay the vortex elongates (weakens) while its 
radial width remains 0{H), so its surface density decreases. 
In addition, for P = 0.1, 0.5 and 1.0 we also observe the ap¬ 
pearance of spiral waves associated with the vortex, which 
may contribute to its dissipation (see below). We will see in 
the next section that this decay after linear growth is very 
different to when the planet potential is kept on. 


5 NON-LINEAR EVOLUTION OF GAP-EDGE 
VORTICES WITH FINITE COOLING TIME 

We now examine long-term simulations of gap-edge vortices 
for P = 0.1, 0.5, 1,5,10. (Additional cases are presented in 
H5.4l when examining vortex lifetimes as a function of p.) The 
planet potential is kept on throughout. We employ a grid 
with {Nr,N^) = (512,1024) in order for these simulations 
to be computationally feasible. We also use a larger disc 
with Tout = 45rin to minimise boundary effects on vortex 
evolution, and apply open boundaries at r = Tin, rout- 

We comment that lower-resolution simulations with 
{Nr,N^) = (256,512) show similar behaviour and trends 
as the high-resolution runs reported below. 



t/Po 


Figure 4. Long term simulations without the planet potential 
after the gap is set up. The m = 1 surface density component, 
non-dimensionalized by the initial m = 0 component, at the outer 
gap edge is shown as a function of the cooling timescale. 



Figure 5. Evolution of the non-dimesnionalized m = 1 surface 
density component for long term simulations with the planet po¬ 
tential kept on. Notice there is vortex growth after the initial 
linear growth, which contrasts to Fig. U where vortices decay in 
the absence of continuous disc-planet interaction. 

5.1 Generic evolution 

The linear growth of the RWI and vortex-formation is fol¬ 
lowed by vortex merging. We now find merging timescales 
independent of P, and by 60Po only one vortex remains. 
The evolution of the amplitude of the m = 1 surface den¬ 
sity component, averaged over r — VpG [2, lOjr^, is shown in 
Fig. [5]for different p. The initial, post-merger vortex ampli¬ 
tude is found to be weaker for longer cooling rates (which 
have smaller linear growth rates). 

In all cases the system remains in a quasi-steady state 
for > SOOPo with a single vortex circulating the outer gap 
edge at the local Keplerian frequency. Fig.[^shows a typical 
plot of the relative surface density perturbation in this state. 
During this stage, the vortex intensifies. This is better shown 
in Fig.[7]as the evolution of Rossby numbers measured at the 
vortex centres. The Rossby number increases in magnitude 
during quasi-steady state, but the maximum |Po| is similar 
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Figure 8 . Average value of the relative surface density pertur¬ 
bation at vortex centres for various cooling times. Initial vortex 
over-densities decrease with cooling rate while vortex growth rates 
increase with cooling rate. 


Figure 6 . Relative surface density perturbation for the /3 = 0.1 
case during quasi-steady state with a single vortex at the outer 
gap edge. The plot for other values of the cooling time $ are 
similar. The decrease in the surface density near r 40rin arises 
from mass loss due to the open boundary condition imposed at 

^out- 



t/Po 


Figure 7. Evolution of Rossby numbers at the centres of the 
vortices formed in discs with different cooling rates. A negative 
Rossby number implies anti-cyclonic motion. Boxcar averaging 
was used to remove contributions from the planet-induced spiral 
shock. 

for all j3: the vortex reaches a characteristic value of Ro « 
—0.35 for ^ = 0.1 and Ro ~ —0.45 for P = 10. 

We find the vortices become significantly over-dense. 
Fig. |S] plots the surface density perturbation measured at 
the vortex centres, showing AS/Eq > 7 for all cases of P 
in qnasi-steady state, and max(AE/So) ~ 11 for /? = 5. 
The maximnm over-density typically increases with longer 
cooling times, despite the vortices are initially weaker at for¬ 
mation with increasing p. The large increase in the surface 
density is due to vortex growth as there is continuous gen¬ 


eration of vorticity by planet-disc interaction. This is sup¬ 
ported by the observation that in the previous simulations 
without the planet, the amplitude of the post-merger vortex 
does not grow (Fig. |3]|. 

Fig.[5]shows that the duration of the quasi-steady state 
varies with the cooling rate: for P = 0.1 and 10, the vortex 
amplitude begins to decay around t SOOPo, while for P = 
1, 5 the decay begins at t ~ 1200Po. This non-monotonic 
dependence suggests that there exists an optimal cooling 
rate to maximise the vortex lifetime. We will discuss this 
issue further in i )5.4l Notice the decay timescale can be long 
with rapid cooling: for P = 0.1 it takes ~ 400Po whereas for 
/3 = 10 it takes ~ lOOPo for the m = 1 amplitude to decay 
significantly after reaching maximum. 

5.2 Additional analysis on vortex decay 

In this subsection we examine the vortex decay observed in 
our simulations in more detail. Fig. show snapshots of the 
vortex for the case = 1. The plots show the surface density 
perturbation and the surface density gradient during quasi¬ 
steady state (t = 700Po), when the m = 1 amplitude begins 
to decrease {t = ISOOPq) and just after the rapid amplitude 
decay (t = 1510Po). 

In quasi-steady (t = 700Po) the vortex is elongated with 
a vortex aspect ratio « 4 , but becomes more compact ap¬ 
proaching a ratio of 2 during its decay {t = ISlOPo). Notice 
in Fig.[5]the appearance of wakes extending from either side 
of the vortex at t = 1300Po. These wakes correlate with 
large gradients in surface density (bottom panel), and are 
first seen in the later half of the quasi-steady state. We find 
the time at which the vortex begins to decay coincides with 
the emergence of these wakes. 

During the quasi-steady state the vortex orbits at r ~ 
Vp + 6rh- We do not see significant vortex migration at this 
stage, since the vortex is located at a surface de nsity max¬ 
imum JPaardekooper. Lesur fc Papaloizoulboiol l. However, 
simultaneous with the appearance of the wakes, we observe 
the vortex begins to migrate inwards to r rp -|- 5rh- 

During quasi-steady state the average value of the sur- 
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Figure 9. The vortex in the case with /3 = 1 during quasi-steady 
state (left), start of decay (middle), and just after the decay in the 
m = 1 amplitude (right). The surface density perturbation (top) 
and the associated surface density gradient (bottom) are shown. 
Wake-like features corresponding to large density gradients are 
found to originate from the vortices during the late phase of their 
quasi-steady states and into dissipation times. This plot is to be 
considered in conjunction with Fig. [5] 


face density gradient along the wakes is |u;sVE/E| ~ 0.4, 
where Ws ~ 0.1 (code units) is a typical length scale of the 
surface density variation across the wake. Just before the 
m = 1 amplitude begins to decrease, we observe this quan¬ 
tity sharply increases to ~ 0.6, and remains around this 
value until the vortex dies out, at which point the associ¬ 
ated Rossby number begins decreasing to zero. After the 
vortex reaches small amplitudes (1 < AE/E), it migrates 
out to r ^ Vp + 6.5rh. 

We also measured large increases in the Mach number 
near the vortex as the m = 1 surface density amplitude 
reaches maximum and begins to decay. Fig. [10] plots the 
Mach number M = \v — Uvor|/cs, where Uvor corresponds 
to the bulk velocity of the vortex around the disc. Values 
in Fig. [To] have been averaged over a region within 2H of 
the vortex centre. During the quasi-steady state the Mach 
number increases steadily, and for all cases M maximizes 
about lOOPo after the start the m = 1 surface density 
amplitude starts to decay. 

Putting the above observations together, we suggest 
that vortex decay (in the m = 1 surface density ampli¬ 
tude) is due to shock formation by the vortex. When the 
vortex reaches large amplitude, it begins to induce shocks 
in the surrounding fluid, as supported by the increase in 
Mach number and the appearance of wakes with large sur¬ 
face density gradients. The vortex may lose energy through 
shock dissipation. In addition, a strong vortex (or shock for¬ 
mation) can smooth out the gap structure that originally 
gave rise to the RWI, which would oppose vortex growth. 
We examine this below. 


5.3 Vortex decay and gap structure 

We find vortex decay modifies the gap structure. Fig. [11] 
shows the gap profile before and after vortex decay for the 



Figure 10. Mach number relative to the vortex, averaged over a 
region within 2H of the vortex centre with respect to time. This 
plot can be compared to the evolution of the vortex amplitude 
shown in Fig. [5] 


^= 1.0 



('■-rp)/r, 


Figure 11. Azimuthally averaged profiles of the relative surface 
density perturbations for = 1.0 before (solid) and after vor¬ 
tex decay (dashed). Alongside drastic reduction in the outer gap 
maxima, vortex decay smooth out the sharp outer gap edge. 


case 13 = 1. The vortex resides around the local surface 
density maximum at the outer gap edge (r ~ rp + 6r^). We 
see that after its amplitude has decayed {t ~ ISOOPo, Fig. 
0, this local surface density maximum is also smoothed out. 

We characterize the smoothness of the outer gap edge 
with a dimensionless gap edge gradient parameter 


JE(t) = 


dr 


(E(t = 0,r)),i, 


(17) 


where Ar = r € [vp, rp-\-Qrh\ is the radial range of averaging, 
spanning from centre of the gap to the radius of the surface 
density maximum. A larger JE characterizes a sharper gap 
edge and larger local surface density maxima. 

A plot of the gradient parameter over time for the /3 = 
1.0 case is shown in Fig. 1121 During vortex decay, the outer 
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Figure 12. Non-dimensional measure of the surface density gra¬ 
dient at the outer gap edge /3 = 1.0. During the vortex quasi¬ 
steady state the gap edge is found to have a large gradient and 
sharp peak while vortex dissipation, which occurs at t ^ IO^Pq 
as seen in Fig. [5] works to smooth out the gap edge. 

gap edge is drastically smoothed out, changing from a value 
of 5E = 1.2 during quasi-steady state to 0.4 after dissipation. 

This can be interpreted as the vortex providing a viscos¬ 
ity; and we measnre a typical alpha viscosity a = 0(10“^) 
associated with the vortex. This acts against gap-opening by 
the planet, and smooths out the outer surface density bump, 
so the condition for the RWI becomes less favonrable. In or¬ 
der to re-launch the RWI, the surface density bump should 
reform. However, this is difficult as there is no more material 
in the planet’s vicinity to clear out (Fig. Hill to form a sur¬ 
face density bump outside the gap. This may explain why 
vortices do not reform again (at least within the simulation 
timescale). After fnll decay the aspect ratio at the outer gap 
edge is h 0.05. 

5.4 Vortex lifetimes as a function of cooling rate 

We now examine vortex lifetimes as function of the imposed 
cooling times. For this study, additional simulations with 
/3 = 0.01, 0.25, 0.75, 2.5, 7.5 were also performed. We define 
the vortex lifetime, tufe, as the time at which the over-density 
of the vortex returns to AE/Eq 1 after reaching maximum 
(which is on the order of the initial over-density associated 
with the gap formation). 

We plot tiife with respect to cooling times in Fig. 1131 We 
also plot tdiss: the time elapsed before the vortex to begins to 
dissipate (when the m — 1 surface density amplitude begins 
to decay); and tMach: the time taken for the average Mach 
number around the vortex to maximise. 

For fast cooling rates (/3 < 1), the vortex lifetime is 
maximized for /3 —>■ 0: we find tufe ~ 2100Po for P = 0.01 and 
decreases to tufe ss 950Po for /3 = 0.25. Note that for very 
small /3, there is significant contribution to the overall vortex 
lifetime due to a long decay timescale. For longer cooling 
times (/3 > 1) the vortex lifetimes maximizes at /) = 2.5 
with tiife « 1650Po- 

We comment here that a locally isothermal simnlation, 
where the disc sonnd-speed is kept constant in time, was also 
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Figure 13. Characteristic timescales associated with vortex evo¬ 
lution, as a function of the cooling parameter /3: time in which the 
m = 1 surface density amplitude begins to decay, tdiss (black); 
time at which the over-density at the vortex centre decreases to 
AS/S ~ 1 after reaching maximum, iufe (red); and tjuach is the 
time taken for the average Mach number around the vortex to 
maximise (blue). 

performed for comparison. In this case we did not observe 
significant vortex decay within the simulation timescale (im¬ 
plying tdiss ^ 2OOOP0). We expect the corresponding vortex 
lifetime to exceed that for /? = 0.01. Inclnding an energy 
equation with rapid cooling (or setting 7 close to unity), 
could still lead to discrepancies with a locally isothermal 
disc. This is due to the advection-creation of specific en¬ 
tropy within the planet’s horseshoe region with the former 
case, thereby affecting the generalized vortensity and there¬ 
fore the instability f 114.31) and snbsequent vortex evolntion. 
Nevertheless, a longer vortex lifetime in a locally isothermal 
disc would be consistent with the above trend of increasing 
vortex lifetimes as /3 —>■ 0. 

In the previous section, we observed that vortices began 
to decay when it starts to induce shocks. We thus suggest 
that the time needed for the vortex to grow to sufficient 
amplitude to induce shocks in the surrounding fluid, which 
may be considered as the duration of the qnasi-steady state 
or tdiss, is an important contribution to the overall vortex 
lifetime. We discnss below some competing factors that may 
result in a non-monotonic dependence of tdiss on the cooling 
rate. 

5.4-1 Factors that lengthen vortex lifetimes 

It has been shown that the amplitude at which the RWI sat- 
ura tes increases with the growth r ate of the linear instabil¬ 
ity llMeheut. Lovelace &: Laill2013l ). Our ‘planet-off’ simnla- 
tions yield slower growth rates with increasing cooling times, 
which suggest weaker vortices are formed initially with in¬ 
creasing p. This is consistent with the present simulations: 
at the beginning of the quasi-steady state {t IOOPq) we 
find the over-density at the vortex centre is AE/Eq = 2.5 
for ,5 = 0.1 and AE/Eq = 1.48 for P = 10. 

The growth of the post-merger single vortex is mediated 
by disc-planet interaction. However, gap-opening becomes 
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more difficult in a hotter disc, and we find the generalised 
vortensity profiles are smoother with increasing (3. This op¬ 
poses the RWI. Furthermore, the vortex should reach larger 
amplitudes to induce shocks on account of the increased 
sound-speed. 

These considerations suggest, with increased cooling 
times, it takes longer for the post-merger vortex grow to suf¬ 
ficient amplitude to induce shocks and dissipate. This factor 
contributes to a longer quasi-steady state with increasing [3. 


5.4-2 Factors that shorten vortex lifetimes 

Notice in Fig. [5] and Fig. [S] the vortex growth during the 
quasi-steady state is actually faster for /3 = 10 than for 
(3 = 5. For example, at f ~ 500Po the vortex with (3 = IQ 
has a larger amplitude than for (3 = 5. This is also reflected 
in Fig. 1101 where the Mach number reaches its maximum 
value sooner for (3 = IQ than for (3 = 5. 

This observation is consiste nt with the RWI being 
favoured by higher temperatures dLi et al.ll200d : Ii^l2012l i 
through the perturbations (as opposed to its effect through 
the set up of the gap profile discussed previously), which 
corresponds to longer cooling times in our case. While our 
‘planet-off’ simulations indicate this is unimportant for the 
linear instability, it may have contributed significantly to the 
vortex growth during quasi-steady state at very long cooling 
times (e.g. (3 = 10). This effect shortens the vortex lifetime 
by allowing it to grow faster and induce shocks sooner. 


6 SUMMARY AND DISCUSSION 

In this paper, we have carried out numerical simulations 
of non-isothermal disc-planet interaction. Our simulations 
were customized to examine the effect of a finite cooling 
time on the stability of gaps opened by giant planets to 
the so-called vortex or Rossby wave instability. To do so, 
we included an energy equation with a cooling term that 
restores the disc temperature to its initial profile on a char¬ 
acteristic timescale t^. We studied the evolution of the gap 
stability as a function of t^. This is a natural extension to 
previous studies of on gap stability, which employ locally 
or strictly isothermal equations of s tate. We considered the 
invis cid limit which favors the RWI dLi et al ] |2009l: [Rr et al.l 
|2014 1 and avoids complications from viscous heating other 
that shock heating. However, this means that the vortex life¬ 
times observed in our simulations are likely longer than in 
realistic discs with non-zero physical viscosity. 

We considered two types of numerical experiments. We 
first used disc-planet interaction to self-consistently set up 
gap profiles, which were then perturbed and evolved without 
further the influence of the planet potential. This procedure 
isolates the effect of cooling on gap stability through the 
set up of the initial gap profile. We find that as the cool¬ 
ing time tc is increased, the gaps became more stable, with 
lower growth rates of non-axisymmetric modes and the dom¬ 
inant azimuthal wavenumber also decreases. This is consis¬ 
tent with the notion that increasing tc leads to higher tem¬ 
peratures or equivalently the disc aspect ratio h, which op¬ 
poses gap-opening by the planet. This means that the gaps 
opened by the planet in a disc with longer tc are smoother 
and therefore more stable to the RWI. 


In the second set of calculations, we included the planet 
potential throughout the simulations and examined the 
long-term evolution of the gap-edge vortex that develops 
from the RWI. The vortex reaches a quasi-steady state last¬ 
ing 0(10®) orbits. Unlike the ‘planet-off’ simulations, in 
which vortices decay after linear growth and merging, we 
find that with the planet potential kept on, the vortex am¬ 
plitude grows during this quasi-steady state, during which 
no vortex migration is observed, until it begins to induce 
shocks, after which the vortex amplitude begins to decay. 

For our main simulations with (3 ^ 0.1, the duration 
of the quasi-steady state increases with increasing cooling 
timescales until a critical value, beyond which this quasi¬ 
steady state shortens again. We find the timescale for the 
vortex to decay after reaching maximum amplitude can be 
long for small (3, which contributes to a long overall vortex 
lifetime with rapid cooling. We do observe vortex migration 
during its decay, which may influence this decay timescale. 

We suggest a non-monotonic dependence of the quasi¬ 
steady state on the cooling timescale (3 can be attributed to 
the time required for the vortex to grow to sufficient am¬ 
plitude to induce shocks in the surrounding fluid, thereby 
losing energy and also smooth out the gap edge. 

For short cooling timescales, the planet is able to open a 
deeper gap which favours the RWI, leading to stronger vor¬ 
tices. For long cooling timescales, we find the vortex grows 
faster during the quasi-stead y state. In acco rdance with pre¬ 
vious stability calculations ( Li et al]l2000l l. we suggest the 
latter is due to the RWI being favoured with increasing disc 
temperature, and that this effect overcomes weaker gap¬ 
opening for sufficiently long cooling times. These competing 
factors imply for both short and long cooling timescales, the 
vortex reaches its maximum amplitude, shock, and begins 
to decay, sooner than intermediate cooling timescales. 

(However, for very rapid cooling, e.g. (3 = 0.01, the 
quasi-steady state is also quite long. This suggests that the 
above effects themselves do not have a simple dependence on 
the cooling timescale when considering /3 —^ 0 and/or that 
other factors become important in this limit. This should be 
investigated in future works.) 

We remark that a non-monoto nic depen d ence of the 
vortex lifetime was also reported by IFu et all (120141 ') , who 
performed locally isothermal disc-planet simulations with 
different values of the disc aspect ratio. In their simulations 
the optimum aspect ratio is h = 0.06. In our simulations, h 
is a dynamical variable, but by analyzing the region where 
the vortex is located {r — Vp € [2,10]r;,), we find for a dimen¬ 
sionless cooing timescale of (3 = 2.5, which has the longest 
vortex lifetime in the presence of moderate cool ing, that 
h ~ 0 .058 on average. Our result is consistent with iFu et al.l 
(l2014lf . 


6.1 Caveats and outlooks 

There are several outstanding issues that needs to be ad¬ 
dressed in future work: 

Convergence. Although lower resolution simulations 
performed in the early stages of this project gave similar 
results (most importantly, the non-monotonic dependence 
of vortex lifetimes on the cooling timescale), we did find the 
lower resolution typically yield longer vortex lifetimes than 
that reported in this paper. This could be due to weaker 
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RWI with low resolution. It will be necessary to perform 
even higher resolution simulations in order to obtain quan¬ 
titatively converged vortex lifetimes. 

Orbital migration. We have held the planet fixed on a 
circular orbit. However, gap-edge vortices are kno wn to ex¬ 
ert s ignificant, oscillatory torques on the planet (iLi et al.l 
I2OO9I I which can lead to complex orbital migration. This will 
likely affect vortex lifetimes as it may alter the planet-vortex 
separati on, as well as lead i ng to direct vortex-plane t inter¬ 
actions (iLin fc Papaloiz^ l2010l : lAtaiee et al. 1 12014) . Thus, 
future simulations should allow the planet to freely migrate. 
Similarly, the role of vortex migration on its lifetime should 
be clarified. 

Cooling model. Our prescription for the disc heat¬ 
ing/cooling is convenient to probe the full range of ther¬ 
modynamic response of the disc. However, in order to cal¬ 
culate vortex lifetimes in actual protoplanetary discs, an 
improved thermodynamics treatment, e.g. radiative cooling 
based on realistic disc temperature, density, opacity models 
etc., should be used in future work. 

Self-gravity. We have ignored dis c self-gravity in this 
study . Based on linear calculations, iLovelace fc Hohlfeldl 
(l2013l) concluded self-gravity to be important for the RWI 
when the Toomre parameter Q < 0{l/h), or Q < 20 for 
h ~ 0.05, as was typically considered in this work. This sug¬ 
gests that self-gravity may affect vortex lifetimes even when 
Q is not small. In particular, given that we observe vortices 
can reach significant over-densities (up to almost an order of 
magnitude), it will be important to include disc self-gravity 
in the future. 

Three-dimensional (3D) effects. A vortex in a 3D disc 
may be subject to secondar y instabilities that destroy them 
dLesur fc Papaloizoul l2009l : iRailton fc Panaloizoul l20l4) . 
This may be an important factor in determining gap-edge 
vortex lifetimes in realistic discs. For example, if these sec¬ 
ondary instabilities sets in before the vortex grows to suffi¬ 
cient amplitude to shock, then the dependence of the vortex 
lifetime on the cooling timescale will be its effect through 
the 3D instability (as opposed to the effect on the RWI it¬ 
self, which is a 2D instability). This problem needs to be 
clarified with full 3D disc-planet simulations. 
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